home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 4 / Amiga Tools 4.iso / grafix / raytracing / raylab / source / intersct.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-01-04  |  12.5 KB  |  422 lines
  1. /*
  2.     name:    intersct.c
  3.  
  4.     Intersection-routines
  5.     ---------------------
  6.  
  7.     These routines will calculate the t-parameter of a line,
  8.     which corresponds to the point where the given line intersects with
  9.     the given primitive. The routines return zero (0) if no intersection
  10.     occured.
  11.  
  12. */
  13.  
  14. #include  <stdio.h>
  15.  
  16. #include  "defs.h"
  17. #include  "extern.h"
  18.  
  19.  
  20. /**********************************************************************
  21.  *
  22.  *    Get (possible) intersection with given object
  23.  *
  24.  **********************************************************************/
  25.  
  26.  
  27. double Intersect_Object(LINE *Line, OBJECT *Object)
  28. {
  29.     LINE    TransformedLine;
  30.     double    t;
  31.     long    i;
  32.     VECTOR    tempv;
  33.  
  34.     CopyPoint(&TransformedLine.Origin,&Line->Origin);
  35.     CopyVector(&TransformedLine.Direction,&Line->Direction);
  36.  
  37.     if(Object->Transform.NumTransforms>0) {
  38.         for(i=Object->Transform.NumTransforms-1L;i>=0L;i--) {
  39.         switch(Object->Transform.Entry[i].Type) {
  40.             case TRANSFORM_SCALE:
  41.             TransformedLine.Origin.x=TransformedLine.Origin.x/Object->Transform.Entry[i].Values.x;
  42.             TransformedLine.Origin.y=TransformedLine.Origin.y/Object->Transform.Entry[i].Values.y;
  43.             TransformedLine.Origin.z=TransformedLine.Origin.z/Object->Transform.Entry[i].Values.z;
  44.             TransformedLine.Direction.x=TransformedLine.Direction.x/Object->Transform.Entry[i].Values.x;
  45.             TransformedLine.Direction.y=TransformedLine.Direction.y/Object->Transform.Entry[i].Values.y;
  46.             TransformedLine.Direction.z=TransformedLine.Direction.z/Object->Transform.Entry[i].Values.z;
  47.             break;
  48.             case TRANSFORM_MOVE:
  49.             TransformedLine.Origin.x=TransformedLine.Origin.x-Object->Transform.Entry[i].Values.x;
  50.             TransformedLine.Origin.y=TransformedLine.Origin.y-Object->Transform.Entry[i].Values.y;
  51.             TransformedLine.Origin.z=TransformedLine.Origin.z-Object->Transform.Entry[i].Values.z;
  52.             break;
  53.             case TRANSFORM_ROTATE:
  54.             NegVector(&tempv,&Object->Transform.Entry[i].Values);
  55.             RotatePoint(&TransformedLine.Origin,&tempv);
  56.             RotateVector(&TransformedLine.Direction,&tempv);
  57.             break;
  58.             case TRANSFORM_NONE:
  59.             default:
  60.             break;
  61.         }
  62.         }
  63.     }
  64.  
  65.     switch(Object->ShapeType) {
  66.         case SHAPE_PLANE:
  67.         t=Intersect_Plane(&TransformedLine, (PLANE *)(Object->Shape));
  68.         break;
  69.         case SHAPE_SPHERE:
  70.         t=Intersect_Sphere(&TransformedLine, (SPHERE *)(Object->Shape));
  71.         break;
  72.         case SHAPE_ELLIPSOID:
  73.         t=Intersect_Ellipsoid(&TransformedLine, (ELLIPSOID *)(Object->Shape));
  74.         break;
  75.         case SHAPE_TRIANGLE:
  76.         t=Intersect_Triangle(&TransformedLine, (TRIANGLE *)(Object->Shape));
  77.         break;
  78.         case SHAPE_BOX:
  79.         t=Intersect_Box(&TransformedLine, (BOX *)(Object->Shape));
  80.         break;
  81.         case SHAPE_DISC:
  82.         t=Intersect_Disc(&TransformedLine, (DISC *)(Object->Shape));
  83.         break;
  84.         case SHAPE_CYLINDER:
  85.         t=Intersect_Cylinder(&TransformedLine, (CYLINDER *)(Object->Shape));
  86.         break;
  87.         default:
  88.         t=0.0;
  89.     }
  90.  
  91.     return(t);
  92. }
  93.  
  94.  
  95.  
  96. /**********************************************************************
  97.  *
  98.  *    These are the actual intersection calculation routines
  99.  *
  100.  **********************************************************************/
  101.  
  102.  
  103. double Intersect_Plane(LINE *Line, PLANE *Plane)
  104. {
  105.     register double    t,pnx,pny,pnz,a,tmp;
  106.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  107.  
  108.     PlaneIntersectionAttempts++;
  109.  
  110.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  111.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  112.     pnx=Plane->Normal.x; pny=Plane->Normal.y; pnz=Plane->Normal.z; a=Plane->a;
  113.  
  114.     tmp=pnx*lvx+pny*lvy+pnz*lvz;
  115.     if(tmp!=0.0) {
  116.         t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
  117.         if(t<EPSILON) { t=0.0; }
  118.         else { PlaneIntersections++; }
  119.     }
  120.     else { t=0.0; }
  121.  
  122.     return(t);
  123. }
  124.  
  125.  
  126.  
  127. double Intersect_Sphere(LINE *Line, SPHERE *Sphere)
  128. {
  129.     double    vx,vy,vz,lx0,ly0,lz0,sx0,sy0,sz0,r;
  130.  
  131.     register double    s,sroot,vxdx,vydy,vzdz,vxsqr,vysqr,vzsqr;
  132.     double    t;
  133.     double    dx,dy,dz;
  134.     double    dxsqr,dysqr,dzsqr;
  135.     double    vsum,t1,t2;
  136.  
  137.     SphereIntersectionAttempts++;
  138.  
  139.     vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
  140.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  141.     sx0=Sphere->Centre.x; sy0=Sphere->Centre.y; sz0=Sphere->Centre.z;
  142.     r=Sphere->r;
  143.  
  144.                         /* These are used more than once, */
  145.     dx=lx0-sx0; dy=ly0-sy0; dz=lz0-sz0;    /* thus much time is saved by not */
  146.     vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;    /* performing the same operation  */
  147.     dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;    /* twice (or even more).          */
  148.     vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
  149.  
  150.     vsum=(vxsqr+vysqr+vzsqr);
  151.  
  152.     s=-(vxdx+vydy+vzdz);
  153.  
  154.     sroot=2*(vxdx*(vydy+vzdz)+vydy*vzdz)+r*r*vsum;
  155.     sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);
  156.  
  157.     if(sroot>=0) {
  158.         sroot=sqrt(sroot);
  159.  
  160.         t1=(s+sroot)/vsum;
  161.         t2=(s-sroot)/vsum;
  162.  
  163.         if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
  164.             t=t1;
  165.             SphereIntersections++;
  166.         }
  167.         else {
  168.             if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
  169.                 t=t2;
  170.                 SphereIntersections++;
  171.             }
  172.             else { t=(double) 0; }
  173.         }
  174.     }
  175.     else t=(double) 0;
  176.  
  177.     return(t);
  178. }
  179.  
  180.  
  181.  
  182. double Intersect_Ellipsoid(LINE *Line, ELLIPSOID *Ellipsoid)
  183. {
  184.     double    vx,vy,vz,lx0,ly0,lz0,ex0,ey0,ez0,a,b,c;
  185.  
  186.     register double    sroot,asqr,bsqr,csqr,vxdx,vydy,vzdz;
  187.     double    t;
  188.     double    dx,dy,dz;
  189.     double    vxsqr,vysqr,vzsqr,dxsqr,dysqr,dzsqr;
  190.     double    r,s,t1,t2;
  191.  
  192.     EllipsoidIntersectionAttempts++;
  193.  
  194.     vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
  195.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  196.     ex0=Ellipsoid->Centre.x; ey0=Ellipsoid->Centre.y; ez0=Ellipsoid->Centre.z;
  197.     a=Ellipsoid->Radius.x; b=Ellipsoid->Radius.y; c=Ellipsoid->Radius.z;
  198.  
  199.     asqr=a*a; bsqr=b*b; csqr=c*c;        /* These are used more than once, */
  200.     dx=lx0-ex0; dy=ly0-ey0; dz=lz0-ez0;    /* thus much time is saved by not */
  201.     vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;    /* performing the same operation  */
  202.     dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;    /* twice (or even more).          */
  203.     vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
  204.  
  205.     s=-((vxdx/asqr)+(vydy/bsqr)+(vzdz/csqr));
  206.  
  207.     sroot=csqr*(vxsqr*bsqr+vysqr*asqr)+vzsqr*asqr*bsqr;
  208.     sroot+=2*(vxdx*(csqr*vydy+bsqr*vzdz)+asqr*vydy*vzdz);
  209.     sroot-=vxsqr*(csqr*dysqr+bsqr*dzsqr);
  210.     sroot-=vysqr*(csqr*dxsqr+asqr*dzsqr);
  211.     sroot-=vzsqr*(bsqr*dxsqr+asqr*dysqr);
  212.  
  213.     if(sroot>=0) {
  214.         sroot=sqrt(sroot)/(a*b*c);
  215.  
  216.         r=((vxsqr/asqr)+(vysqr/bsqr)+(vzsqr/csqr));
  217.  
  218.         t1=(s+sroot)/r;
  219.         t2=(s-sroot)/r;
  220.  
  221.         if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
  222.             t=t1;
  223.             EllipsoidIntersections++;
  224.         }
  225.         else {
  226.             if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
  227.                 t=t2;
  228.                 EllipsoidIntersections++;
  229.             }
  230.             else { t=(double) 0; }
  231.         }
  232.     }
  233.     else t=0.0;
  234.  
  235.     return(t);
  236. }
  237.  
  238.  
  239. double Intersect_Triangle(LINE *Line, TRIANGLE *Triangle)
  240. {
  241.     double    t=0.0,ttmp,pnx,pny,pnz,a;
  242.     double    lvx,lvy,lvz,lx0,ly0,lz0,angle1,angle2,angle3;
  243.     VECTOR    v1,v2,v3;
  244.     POINT    ip;
  245.  
  246.     TriangleIntersectionAttempts++;
  247.  
  248.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  249.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  250.     pnx=Triangle->Plane.Normal.x; pny=Triangle->Plane.Normal.y; pnz=Triangle->Plane.Normal.z; a=Triangle->Plane.a;
  251.  
  252.     ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/(pnx*lvx+pny*lvy+pnz*lvz);
  253.     if(ttmp>EPSILON) {
  254.         ip.x=Line->Origin.x+ttmp*Line->Direction.x;    /* ip = intersection point */
  255.         ip.y=Line->Origin.y+ttmp*Line->Direction.y;
  256.         ip.z=Line->Origin.z+ttmp*Line->Direction.z;
  257.  
  258.         if((ip.x>Triangle->Min.x)&&(ip.x<Triangle->Max.x)&&(ip.y>Triangle->Min.y)&&(ip.y<Triangle->Max.y)&&(ip.z>Triangle->Min.z)&&(ip.z<Triangle->Max.z)) {
  259.         CreateVector(&v1,Triangle->Corners[1].x-Triangle->Corners[0].x,Triangle->Corners[1].y-Triangle->Corners[0].y,Triangle->Corners[1].z-Triangle->Corners[0].z);
  260.         CreateVector(&v2,Triangle->Corners[2].x-Triangle->Corners[0].x,Triangle->Corners[2].y-Triangle->Corners[0].y,Triangle->Corners[2].z-Triangle->Corners[0].z);
  261.         angle1=VectorsAngle(&v1,&v2);
  262.         CreateVector(&v3,ip.x-Triangle->Corners[0].x,ip.y-Triangle->Corners[0].y,ip.z-Triangle->Corners[0].z);
  263.         angle2=VectorsAngle(&v1,&v3);
  264.         angle3=VectorsAngle(&v2,&v3);
  265.         if((angle2+angle3)<=(angle1+EPSILON)) {
  266.             CreateVector(&v1,Triangle->Corners[1].x-Triangle->Corners[2].x,Triangle->Corners[1].y-Triangle->Corners[2].y,Triangle->Corners[1].z-Triangle->Corners[2].z);
  267.             CreateVector(&v2,Triangle->Corners[0].x-Triangle->Corners[2].x,Triangle->Corners[0].y-Triangle->Corners[2].y,Triangle->Corners[0].z-Triangle->Corners[2].z);
  268.             angle1=VectorsAngle(&v1,&v2);
  269.             CreateVector(&v3,ip.x-Triangle->Corners[2].x,ip.y-Triangle->Corners[2].y,ip.z-Triangle->Corners[2].z);
  270.             angle2=VectorsAngle(&v1,&v3);
  271.             angle3=VectorsAngle(&v2,&v3);
  272.             if((angle2+angle3)<=(angle1+EPSILON)) {
  273.                 t=ttmp;
  274.                 TriangleIntersections++;
  275.             }
  276.         }
  277.         }
  278.     }
  279.  
  280.     return(t);
  281. }
  282.  
  283.  
  284. double Intersect_Box(LINE *Line, BOX *Box)
  285. {
  286.     double    t=0.0,told=-1.0,ttmp,pnx,pny,pnz,a;
  287.     double    lvx,lvy,lvz,lx0,ly0,lz0;
  288.     int    i;
  289.     POINT    ip;
  290.  
  291.     BoxIntersectionAttempts++;
  292.  
  293.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  294.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  295.     for(i=0;i<6;i++) {
  296.         pnx=Box->Planes[i].Normal.x; pny=Box->Planes[i].Normal.y; pnz=Box->Planes[i].Normal.z; a=Box->Planes[i].a;
  297.  
  298.         ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/(pnx*lvx+pny*lvy+pnz*lvz);
  299.         if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
  300.             switch(i) {
  301.                 case 0:
  302.                 case 1:
  303.                 ip.y=Line->Origin.y+ttmp*Line->Direction.y;
  304.                 ip.z=Line->Origin.z+ttmp*Line->Direction.z;
  305.                 if((ip.y>=Box->Corners[0].y)&&(ip.y<=Box->Corners[1].y)&&(ip.z>=Box->Corners[0].z)&&(ip.z<=Box->Corners[1].z)) {
  306.                     told=ttmp;
  307.                 }
  308.                 break;
  309.                 case 2:
  310.                 case 3:
  311.                 ip.x=Line->Origin.x+ttmp*Line->Direction.x;
  312.                 ip.z=Line->Origin.z+ttmp*Line->Direction.z;
  313.                 if((ip.x>=Box->Corners[0].x)&&(ip.x<=Box->Corners[1].x)&&(ip.z>=Box->Corners[0].z)&&(ip.z<=Box->Corners[1].z)) {
  314.                     told=ttmp;
  315.                 }
  316.                 break;
  317.                 case 4:
  318.                 case 5:
  319.                 ip.x=Line->Origin.x+ttmp*Line->Direction.x;
  320.                 ip.y=Line->Origin.y+ttmp*Line->Direction.y;
  321.                 if((ip.x>=Box->Corners[0].x)&&(ip.x<=Box->Corners[1].x)&&(ip.y>=Box->Corners[0].y)&&(ip.y<=Box->Corners[1].y)) {
  322.                     told=ttmp;
  323.                 }
  324.                 break;
  325.             }
  326.         }
  327.     }
  328.     if(told>EPSILON) {
  329.         t=told;
  330.         BoxIntersections++;
  331.     }
  332.  
  333.     return(t);
  334. }
  335.  
  336.  
  337. double Intersect_Disc(LINE *Line, DISC *Disc)
  338. {
  339.     register double    t,pnx,pny,pnz,a,tmp;
  340.     double    lvx,lvy,lvz,lx0,ly0,lz0,dx,dy,dz;
  341.  
  342.     DiscIntersectionAttempts++;
  343.  
  344.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  345.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  346.     pnx=Disc->Plane.Normal.x; pny=Disc->Plane.Normal.y; pnz=Disc->Plane.Normal.z; a=Disc->Plane.a;
  347.  
  348.     tmp=pnx*lvx+pny*lvy+pnz*lvz;
  349.     if(tmp!=0.0) {
  350.         t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
  351.         if(t<EPSILON) { t=0.0; }
  352.         else {
  353.             dx=(Line->Origin.x+Line->Direction.x*t)-Disc->Centre.x;
  354.             dy=(Line->Origin.y+Line->Direction.y*t)-Disc->Centre.y;
  355.             dz=(Line->Origin.z+Line->Direction.z*t)-Disc->Centre.z;
  356.             if((dx*dx+dy*dy+dz*dz)>(Disc->r*Disc->r)) { t=0.0; }
  357.             else { DiscIntersections++; }
  358.         }
  359.     }
  360.     else { t=0.0; }
  361.  
  362.     return(t);
  363. }
  364.  
  365.  
  366. double Intersect_Cylinder(LINE *Line, CYLINDER *Cylinder)
  367. {
  368.     register double    t=0.0,t1,t2,t3,t4,rsqr,lx0lvx,ly0lvy;
  369.     double    lvx,lvy,lvz,lx0,ly0,lz0,tmp,tmp2,tmp3;
  370.     POINT    ip1,ip2;
  371.  
  372.     CylinderIntersectionAttempts++;
  373.  
  374.     lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
  375.     lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
  376.     rsqr=Cylinder->r*Cylinder->r;
  377.     
  378.     lx0lvx=lx0*lvx; ly0lvy=ly0*lvy;
  379.     tmp=2.0*lx0lvx*ly0lvy-lvx*lvx*ly0*ly0-lvy*lvy*lx0*lx0+lvx*lvx*rsqr+lvy*lvy*rsqr;
  380.  
  381.     if(tmp>=0.0) {
  382.         tmp=sqrt(tmp);
  383.         tmp2=-(lx0lvx+ly0lvy);
  384.         tmp3=1.0/(lvx*lvx+lvy*lvy);
  385.         t1=(tmp2+tmp)*tmp3;
  386.         t2=(tmp2-tmp)*tmp3;
  387.         ip1.x=lx0+lvx*t1; ip1.y=ly0+lvy*t1; ip1.z=lz0+lvz*t1;
  388.         ip2.x=lx0+lvx*t2; ip2.y=ly0+lvy*t2; ip2.z=lz0+lvz*t2;
  389.  
  390.         tmp=ip1.x*Cylinder->Discs[0].Plane.Normal.x+ip1.y*Cylinder->Discs[0].Plane.Normal.y+ip1.z*Cylinder->Discs[0].Plane.Normal.z;
  391.         if(tmp>Cylinder->Discs[0].Plane.a) { t1=0.0; }
  392.         else {
  393.             tmp=ip1.x*Cylinder->Discs[1].Plane.Normal.x+ip1.y*Cylinder->Discs[1].Plane.Normal.y+ip1.z*Cylinder->Discs[1].Plane.Normal.z;
  394.             if(tmp>Cylinder->Discs[1].Plane.a) { t1=0.0; }
  395.         }
  396.         tmp=ip2.x*Cylinder->Discs[0].Plane.Normal.x+ip2.y*Cylinder->Discs[0].Plane.Normal.y+ip2.z*Cylinder->Discs[0].Plane.Normal.z;
  397.         if(tmp>Cylinder->Discs[0].Plane.a) { t2=0.0; }
  398.         else {
  399.             tmp=ip2.x*Cylinder->Discs[1].Plane.Normal.x+ip2.y*Cylinder->Discs[1].Plane.Normal.y+ip2.z*Cylinder->Discs[1].Plane.Normal.z;
  400.             if(tmp>Cylinder->Discs[1].Plane.a) { t2=0.0; }
  401.         }
  402.  
  403.         t=0.0;
  404.  
  405.         if((t1<EPSILON)||(t2<EPSILON)) {
  406.             t3=Intersect_Disc(Line, &Cylinder->Discs[0]);
  407.             t4=Intersect_Disc(Line, &Cylinder->Discs[1]);
  408.             if(t1>EPSILON) { t=t1; }
  409.             if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) { t=t2; }
  410.             if((t3>EPSILON)&&((t3<t)||(t<EPSILON))) { t=t3; }
  411.             if((t4>EPSILON)&&((t4<t)||(t<EPSILON))) { t=t4; }
  412.         }
  413.         else {
  414.             if(t1>EPSILON) { t=t1; }
  415.             if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) { t=t2; }
  416.         }
  417.         if(t>EPSILON) { CylinderIntersections++; }
  418.     }
  419.  
  420.     return(t);
  421. }
  422.